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Abstract 

New, improved curve fits for the thermodynamic 
properties of equilibrium air have been developed. 
The curve fits are for pressure, speed of sound, tem- 
perature, entropy, enthalpy, density, and internal en- 
ergy. These curve fits can be readily incorporated 
into new or existing computational fluid dynamics 
codes if “real-gas” effects are desired. The curve fits 
are constructed from Grabau-type transition func- 
tions to model the thermodynamic surfaces in a 
piecewise manner. The accuracies and continuity of 
these curve fits are substantially improved over those 
of previous curve fits. These improvements are due 
to the incorporation of a small number of additional 
terms in the approximating polynomials and care- 
ful choices of the transition functions. The ranges 
of validity of the new curve fits are temperatures up 
to 25000 K and densities from 10 -7 to 10’^ amagats. 

Introduction 

Under subsonic flight conditions, air may be 
treated as an ideal gas composed of rigid rotating 
diatomic molecules. The thermodynamic properties 
of such a gas are well known. However, under hyper- 
sonic flight conditions, air may be raised to tempera- 
tures at which the molecules can no longer be treated 
as rigid rotators. Thus, there is a very real need for 
the thermodynamic and transport properties of equi- 
librium air for the computation of flow fields around 
bodies in high-speed flight. The references discussed 
below are representative of the various approaches 
for obtaining thermodynamic properties, but the list 
is by no means complete. 

The thermodynamic properties of equilibrium 
air were calculated with good confidence as early 
as 1950. The earliest approach to compiling these 
properties was to present the information in the form 
of tables or charts (refs. 1 to 4). 

Subsequently, equilibrium air thermodynamic 
properties became available in the form of FOR- 
TRAN computer programs. These programs can be 
broadly divided into two categories. The first cat- 
egory consists of programs that compute the equi- 
librium composition and thermodynamic properties 
using a harmonic-oscillator rigid-rotator model for 
the various component species of the gaseous mix- 
ture. Programs (refs. 5 to 8) were developed for the 
calculation of equilibrium properties of specific gas 
mixtures or of arbitrary chemical systems. 

The second category of computer programs, which 
includes the present work, consists of programs that 
determine the thermodynamic properties of equilib- 
rium air in a noniterative fashion using either in- 
terpolation or polynomial approximation techniques 


(refs. 9 to 16). Typically, the sources of data for these 
programs are references 1 to 4. One such program, 
NASA RGAS (based on ref. 5), was an improvement 
over other sources of thermodynamic properties in 
terms of accuracy and range of validity. For this rea- 
son it is still widely used. The major shortcoming of 
the RGAS program is that the table lookup of coef- 
ficients for the cubic interpolation makes it too cum- 
bersome and time-consuming to be efficiently used 
on an advanced computer. 

Tannehill and associates (refs. 10, 15, and 16) de- 
veloped simplified curve fits for the thermodynamic 
and transport properties of equilibrium air with the 
same range of validity as the NASA RGAS program. 
The curve fits were constructed through the use of 
Grabau-type transition functions in a manner sim- 
ilar to that of reference 11. In forming these curve 
fits, as many as five Grabau-type transition functions 
were joined with the perfect-gas equation of state. 

One of the major shortcomings of the curve fits 
of references 10, 15, and 16 is the lack of continuity 
across the boundaries between the transition func- 
tions. As a consequence, numerical difficulties were 
sometimes encountered when these curve fits were 
employed in iterative flow-field computations. The 
primary objective of the present research was to al- 
leviate this difficulty. At the same time, an attempt 
was made to improve the accuracy of the curve fits 
through incorporation of a small number of addi- 
tional terms which would not significantly increase 
the computation time. 

Through careful choice of the Grabau-type tran- 
sition functions and use of complete bicubic polyno- 
mials, curve fits for pressure, speed of sound, temper- 
ature, entropy, enthalpy, density, and internal energy 
were developed and are presented herein. These 
curve fits are based on the NASA RGAS data and 
have the same ranges of validity, namely, tempera- 
tures up to 25 000 K and densities from 10 -7 to 10 3 
amagats (p/p Q ). 

Symbols 

a speed of sound, m/s 

e specific internal energy, m 2 /s 2 

h specific enthalpy, m 2 /s 2 

p pressure, N/m 2 

R gas constant, 287.06 m 2 /s 2 -K 
s specific entropy, m 2 /s 2 -K 

T temperature, K 

7 = h/e 



p density, kg/nr* 

Subscript: 

o reference conditions at 1 at in and 

273.15 K 

Behavior of Air at High Temperature 

When a gas composed of polyatomic molecules is 
heated to high temperatures, its composition changes 
as a result of the chemical reactions which take place. 
Such a situation exists behind the shock wave which 
envelops a vehicle entering the atmosphere of the 
Earth. As a result of the change in chemical com- 
position. the thermodynamic properties of the gas 
also change. When the temperature of the gas is 
raised appreciably higher than -the temperature at 
which dissociation reactions begin to occur, the elec- 
trons receive energy quanta because of the collisions 
between atoms. If the temperature, and hence the 
kinetic energy of the atoms, is high enough so that 
electrons are removed from their orbits, ionization of 
the gas takes place. The effects of dissociation and 
ionization of the gas on its thermodynamic properties 
are often referred to as “real-gas*' effects. 

At room temperature, the volumetric composition 
of air is about 78 percent diatomic nitrogen, 21 per- 
cent diatomic oxygen, and about 1 percent argon 
and traces of carbon dioxide. When the temperature 
of air is raised above room temperature, deviations 
from perfect- gas behavior occur; that is, the vibra- 
tional mode of the molecules becomes excited, disso- 
ciation of both oxygen and nitrogen molecules occurs 
(although at different temperatures), nitric oxide is 
formed, and so forth. The chemical composition of 
air for densities lying between 10" 2 and 10 times nor- 
mal air density is approximately divisible into the 
following regimes: 

1 T < 2500 K. The chemical composition is sub- 
stantially that at room temperature. 

2. 2500 < T < 4000 K. This is the oxygen disso- 
ciation regime; no significant nitrogen dissocia- 
tion occurs; some NO is formed. 

3. 4000 < T < 8000 K. This is the nitrogen dis- 
sociation regime; oxygen fully dissociates. 

4. T > 8000 K. Ionization of the atomic 
constituents occurs. 

Sources of Equilibrium Air Properties 

Tin' following discussion is intended to summarize 
the available mechanisms for determining equilib- 
rium air properties. The cited references are not in- 
tended as a complete compilation but serve only as a 
list ing typical of the various methods for determining 
the properties. 
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Prior to l%(). methods for determining equilib- 
rium air properties were available only in summary 
form as tables or charts. The sources for information 
were the calculations of Gilmore (ref. 1), Hilsenrath 
and Beckett (ref. 2). Hansen (ref. 3), and Moeckel 
and Weston (ref. 4). In reference 3, data for com- 
pressibility factor, enthalpy, speed of sound, specific 
heat, Prandtl number, and the coefficients of viscos- 
ity and conductivity are presented as functions of 
temperature and pressure. 

Eventually, the calculation of equilibrium air 
properties was possible through the use of FOR- 
TRAN computer programs, which can be divided 
broadly into two categories. The first category 
consists of programs that compute the equilibrium 
composition and thermodynamic properties using a 
harmonic-oscillator rigid-rotator model for the vari- 
ous component species of the gaseous mixture. Bai- 
ley (ref. 5) developed computer programs which 
used the temperature, density, and molar concen- 
trations of the various constituent species to calcu- 
late the pressure, gas constant, enthalpy, entropy, 
specific heats, and coefficient of thermal conductiv- 
ity. These properties were computed for a 9-species 
model as well as an 1 1 -species model of equilibrium 
air. Zeleznik and Gordon (ref. 0) developed a so- 
phisticated computer program, improved later by 
Gordon and McBride (ref. 7), which computed the 
chemical equilibrium composition of complex chemi- 
cal systems given the constituent species and one of 
five possible pairs of thermodynamic state combina- 
tions. Also, a 27- react ion equilibrium air program 
was developed by Miner et al. (ref. 8). 

The second category of computer programs con- 
sists of programs that determine the thermodynamic 
propert ies of equilibrium air in a noniterative fashion 
using either of 1 interpolation-of-polynornial ap- 
proximation techniques. Lomax and Inouye (ref. 9) 
developed FORTRAN programs to determine the 
speed of sound, enthalpy, temperature, and entropy 
as functions of either pressure and density or pres- 
sure and entropy. Their programs used a 9-point 
spline interpolation and required a lookup of over 
10 000 tabulated values. The programs developed at 
NASA Ames Research Center in references 5 and 9 
eventually evolved into the NASA RGAS program. 
The NASA RGAS program employs a cubic inter- 
polation technique, with the associated table lookup 
of cubic coefficients, to compute the enthalpy, tem- 
perature, entropy, and speed of sound of 13 different 
gas mixtures, including equilibrium air as functions 
of either pressure and density, or pressure and en- 
tropy. The NASA RGAS program was modified by 
Tannehill and Molding (ref. 10) to allow internal en- 
ergy and density to be used as independent variables 
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for “time-dependent” flow calculations. The major 
shortcoming of the RGAS program is that the ta- 
ble lookup of coefficients for the cubic interpolation 
makes it too cumbersome and time-consuming to be 
efficiently employed on an advanced computer. 

Among the first to develop programs which 
approximated the thermodynamic properties 
as self-contained closed-form expressions was 
Grabau (ref. 11). He outlined a systematic tech- 
nique of modeling the thermodynamic properties 
with polynomial expressions containing exponential 
t ransitions. Using this technique, he determined the 
enthalpy, entropy, speed of sound, and compressibil- 
ity of equilibrium air as functions of pressure and 
density in the form of closed-form expressions (curve 
fits). Using Grabau’s technique, Lewis and Burgess 
(ref. 12) obtained empirical equations for the density, 
enthalpy, speed of sound, and compressibility factor 
of air as functions of pressure and entropy. How- 
ever, these curve fits had a range of validity only up 
to 15 000 K and a pressure range of 0.1 to 1.0 atm. 
The method of reference 1 1 was also employed by 
Barnwell (ref. 15 ) to curve fit 7 as a function of in- 
ternal energy and density and temperature as a func- 
tion of pressure and density for equilibrium air. Vie- 
gas and Howe (ref. 14 ) developed programs for the 
density, temperature, viscosity, and Prandtl number 
of equilibrium air as functions of pressure and en- 
thalpy in the form of curve fits using least squares 
and Ohebyshev polynomial fitting. Tannehill and 
associates (refs. 10 , 15 , and 16 ) developed simpli- 
fied curve fits for the thermodynamic and transport 
properties of equilibrium air with the same range of 
validity as the NASA RGAS program. These curve 
fits included pressure, temperature, speed of sound, 
and coefficients of viscosity and thermal conductivity 
as functions of internal energy and density; also in- 
cluded were temperature and enthalpy as functions of 
pressure and density. The curve fits were constructed 
using Grabau-type transition functions in a manner 
similar to that of reference 1 1 . In forming these curve 
fits, as many as five Grabau-type transition functions 
were joined with the perfect-gas equation of state. 

Construction of Curve Fits 

Typical Curve Forms 

In the flow calculations of air in thermody- 
namic equilibrium, it becomes important to know 
the various thermodynamic properties as functions 
of a pair of independent state variables. In or- 
der to illustrate the spatial behavior of these ther- 
modynamic surfaces, a typical curve is examined 
here in some detail. The nature of the thermo- 
dynamic surface, with the plausible reasons for its 


undulating behavior, provides a qualitative insight 
into the choice of the approximating functions. Fig- 
ure 1 shows the function 7 plotted with respect to 
loglo(pM;) ~ logio (p/po) at, a density of 1 CT 7 ama- 
gats. Also shown are the various segments into which 
the curve may be divided, as indicated by A, AA, B, 
C, and D. These segments are basically quadratic or 
linear curves which are joined together by transition 
curves. Two types of transition curves appear in fig- 
ure 1, and these are illustrated in figures 2 and 3. 
Figure 2 shows a transition function which passes 
through a point of inflection and is referred to as a 
transition with inflection. Figure 3 illustrates the sec- 
ond type of transition, which is one without a point 
of inflection. Figure 1 shows that 7 goes through 
three distinct transitions with inflections. According 
to reference 3 . there is a definite correlation between 
these three transitions and the change in chemical 
composition of the air as the temperature increases: 
the first transition, from AA to B, is due to the oxy- 
gen dissociation reaction; the second, from B to C, 
is due to the nitrogen dissociation; and the third, 
from C to D, is due to the ionization reactions. 

In addition to the three transitions with inflec- 
tions in figure 1, there appears to be a relatively in- 
significant transition without, an inflection between 
curves A and A A. Also, after a careful examination 
of segment D, it appears that it may actually be part 
of an incomplete transition with a point of inflection. 

The term 7 is plotted as a function of logio(p/Po)~ 
log 1() {p/po) for various densities in figure 4, which 
includes the curve fit of figure 1. As the density 
increases, pieces of the curve near C and D disap- 
pear until only a part of the transition into C re- 
mains at 10 '* amagats. The reason for this is that the 
compressibility factor decreases steadily as the den- 
sity is increased isothermally. Hence, it, also follows 
that, isothermal points move rapidly along the curve 
from D toward C as the density increases. Figure 4 
provides an idea of the complexity of the problem of 
devising a practical method of modelling the collapse 
of the lower segments with increasing density. There 
appears to be a tendency for transitions with inflec- 
tions to convert to transitions without inflections as 
the density increases. Reference 1 suggests that this 
conversion might be correlated with the simultane- 
ous, abrupt increases of the concentrations of ionized 
oxygen and nitrogen atoms and of ionized nitrogen 
molecules. 

As a consequence of the above discussion, one 
is motivated to model the thermodynamic surface, 
in a piecewise manner, with biquadratic or bicubic 
polynomials joined together by exponential transi- 
tion functions with or without points of inflection. 
This is the procedure adopted in the present study. 
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I'igure 2. Transition curve with inflection. 


Figure .’5. Transition curve without inflection. 
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Transition Regions 

The basic forms of the variables 7 and logio(T / T () ), 
plotted at constant densities as functions of 
1 o Kio(p/p») “ log 10 (P/Po), are shown in figures 4 
and 5 . As mentioned previously, these curves ex- 
hibit segments of linear or quadratic functions suc- 
cessively connected by transition functions, which 
are asymptotic at both ends, and may or may not 
include points of inflection. The fact that at least 
some of these transitions can be attributed to dis- 
sociation phenomena suggests the use of exponential 
distribution functions. 

Following the method outlined by Grabau (ref. 11), 
one has a choice of two kernel transition functions. 
The first is the Fermi-Dirac function 


1 + exp(Arz) 


which gradually diminish as the density increases and 
then change to even transitions. There are two ways 
of applying each of these transition functions. When 
the path of a curve appears to move from one straight 
line to another, there is an offset present which can 
be calculated in the direction of either of the vari- 
ables. For accuracy it appears to be better to view 
the transition in terms of the smaller offset. Both 
ways of viewing the offsets involve the choice of a 
baseline. The use of the offset in the y-direction sim- 
plifies this choice since the z-axis serves as a natural 
baseline. 

Consider the problem of determining the equation 
of a curve consisting of two linear segments connected 
by an odd transition function (fig. 8). The lower and 
upper line segments are given by 

yi=(i[X + bi ( 3 ) 


which represents a transition between the levels zero 
and unity, where the direction and rate of the transi- 
tion depend on the sign and the numerical magnitude 
of the exponential constant k. The numerator defines 
the upper level of the transition and may take on a 
variety of forms. In figure 6 the upper level of the 
transition is a straight line inclined to the horizontal, 
while the lower level is the z-axis. The transitions in 
figure (i have points of inflection and, in the termi- 
nology of Grabau (ref. 11), are referred to as odd 
transitions. 

The second type of transition function is the 
kernel of the Bose-Einstein distribution function 


1 - exp(kx) v ; 

which provides transitions leading from one function 
to another without a point of inflection and is ob- 
tained by merely changing the sign before the expo- 
nential term in the denominator of the Fermi-Dirac 
function. The transition function given by equa- 
tion (2) is termed an even transition. Figure 7 illus- 
trates two t ransitions of this kind between the z-axis 
and the line y = z, where (as before) the directions 
and rates of the transitions are governed by the sign 
and magnitude of the exponential constant k. It is 
important to note that the expression for an even 
transition becomes an indeterminate form when z is 
equal to the z-coordinate of the point of intersection 
of t he two lines bounding the transition. 

In the current study, each of the thermodynamic 
curves is approximated with quadratic or incomplete 
cubic segments connected by odd and even transi- 
tions as described above. Almost without exception, 
at low densit ies all the curves undergo odd transitions 


and 


V '2 = «2* + &2 


The y offset is their difference: 


2/2 - 2/1 = (°2 - a i)x + (b 2 - b \ ) 


(4) 

(5) 


which becomes the numerator of the transition func- 
tion. The remaining constants of the transition func- 
tion can be found graphically by drawing three lines 
between y\ and y^. The median line is given by 


Vo = 


y i + m 
2 


(6) 


Let y a be the median line between y G and y\ and y i 
be the corresponding median line between y a and y*)- 
The center of the transition, (z fJ ,y 0 ), is the point at 
which the transition crosses the median line y () . The 
desired transition function is then of the form 


_ , («2 - a \)x+ (62 “ M 

2/ — 2/1 H 7—. 777 vj — 

1 + exp[/c(z - x 0 )\ 


(7) 


The exponential constant k is found from the coor- 
dinates x a and z^ at which the transition intersects 
the lines y a and y/ r Specifically, for the intersection 
with the line y a , 


1 1 

1 + exp [k{x a - *o)] 4 


so that 

exp[fc(z„ - x 0 )} = 3 

Solving for k yields 


k = 


In 3 
x a — x 0 


( 8 ) 

( 9 ) 

(10) 
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Figure (). Two odd transition functions. 



Figure 7. Two even transition functions. 



Figure 8. Construction of an odd transition function. 





From the intersection of y with y b we get 


This procedure obviously yields two numerical values 
for the constant k. However, they are substantially 
alike in most instances. 

The determination of the constants of an even 
transition is simpler. In terms of the y offset, such a 
transition can be written in the form 


, _ u(t - ,r () ) 

1 - exp[fc(x - x„)j 

where x () is the ^-coordinate of the point of intersec- 
tion of the two lines bounding the transition. The 
value of the exponential constant k follows from the 
coordinate y Q at .r = x G . Since the expression for y is 
an indeterminate form at this point, its value is given 
bv the ratio of the derivatives of the numerator and 
of the denominator at this point: 


Vo = — n 1 y (13) 

* hm - {l - ex p[Ar(x-T„)]}} 

which gives 

<M) 

This approach for determining the constants 
of the Grabau-type transitions is extended in the 
present work to approximate transitions in two in- 
dependent variables. The kernel of an odd transition 
function in three dimensions is 


1 

1 + exp(a 0 4 a\ x 4 a 2 y 4 a;\zy) 
which is essentially an alternate form of 


1 

1 + exp[k{x - x 0 )(y - y„) 


(15) 


(16) 


In order to ensure an accurate and smooth transition 
from /i(i\ y) to f 2 (x,y), we require the quadratic 
expression (a 0 4 ayx 4 a 2 y 4 a : yxy) to behave as 
follows. At the lower left corner point (x a , y c ) the 
quadratic expression should have a large positive 
value so that f(x.y) ^ fi(z,y). At the upper right 
corner point (x/,, </,/) the quadratic expression should 
have a large negative value in order to ensure that 
/(•G2/) - hUryy)* At the midpoints of the left 
and right boundaries, [x a , (y r + y d )/ 2] and [x b , (y c + 
y ( l)/'2\, respectively, the quadratic expression should 
be zero so that 


/( x y) % hi£iyl±M^y} 

These conditions yield the following four linear equa- 
tions: 


«o 4 a i s <i 4 a 2 y c 4 a:yx a y c - 4 k (18) 

«() + (l i *b + a 2 Vd + a*x b y d = -k (19) 

<*>0 + a \ x a i + a 2 (y r 4 // f /)/24a3X a (y r 4y f /)/2 = 0 (20) 

a 0 4 fl i Xf, 4 a 2 (yr 4 y ( { )/2 4 a^x b (y r 4 y^/)/2 = 0 (21) 

where k is a positive constant (typically 20 < k < 25) 
chosen such that exp(fc) and exp(— k) do not yield 
overflow and underflow conditions on a computer. 
The constants a<) to a can now be obtained in 
a straightforward manner from the system of four 
linear equations in four unknowns (eqs. (18) to (21)). 

The above met hod of obtaining the Grabau-type 
transition functions proved quite accurate in ensur- 
ing a negligible mismatch in the dependent variable 
over the boundaries of adjoining subregions. It is a 
merit of this stepwise method of constructing em- 
pirical equations that any part can be removed for 
corrections without disturbing the surface approxi- 
mation as a whole. 

Equations of the Curve Fits 


Equation (15) is more convenient for determining the 
values of the constants «o to as dictated by the 
behavior imposed on the transition function. The 
general technique of determining the values of these 
constants differs from the approach outlined earlier 
and is as follows. The boundaries of the transition in 
the directions of the two independent variables are 
x (i < x < x b and y c < y < y (i . If f x (x. y) and f 2 (x, y) 
are the two surfaces limiting the transition function 
/(*,«/), then 


MfiV) ~ fi( x , y) 

1 + exp(a 0 + a x x + a 2 y + a^xy) 


(17) 


The curve fits lor the various thermodynamic 
properties are constructed through use of Grabau- 
type transition functions, as described previously. 
The general form of these curve fits can be written 
as 


z{x } y) = }\{x,y) + 


/ 2 M ~ /i(s,y) 

1 4 exp(Ar 0 + k\x + k 2 y 4 k^xy) 


( 22 ) 


where, in general. 


/ 1 (*. y) = Pi + Pi-r + Psy + P4 x y + p a* 2 + pc,y 2 
+ P 7 x ~y + psxy 2 + pgx 3 + p 10 j/ 3 (23) 


f(x,y) = /1 (i, y) + 



and 


computer program. 


/ 2 (x,y) - h( x ,y) = pii + p i2* + p i3y + p\4 x v 
4 P15^ 2 + P16J/ 2 + Pl7^ 2 y 
4- p\$xy 2 + pigx 3 4 p2oy 3 (24) 


p = p(e,p) 

For the correlation of p = p(e,p), the ratio 7 = 
///c is curve fitted as a function of e and p so that p 
can he calculated from 


The coefficients A’o to in the denominator of the 
transit ion function in equation (22) are determined 
by the technique outlined in the preceding section. 
The coefficients p\ to p 2 o in equations (23) and (24) 
are determined by the actual curve fitting of the data 
from the NASA RGAS program. The exact location 
and number of these data points over the curve fit. 
domain determines the accuracy of the curve fits. 
The points are clustered near the boundaries of the 
domain and the middle region of the transition in 
order to ensure continuity at the boundaries and 
accuracy within the domain. The data from the 
NASA RGAS program are fitted to the equations 
of the curve fits by the method of least squares. A 
multiple linear regression technique (ref. 17) is used 
to determine the coefficients p \ to P20* 

The general form of the curve fit for each ther- 
modynamic property is described below. As in ref- 
erences 10 and 15, for each of the curve fits where 
density is one of the independent variables, the range 
of p is subdivided into three separate regions, with 
different coefficients being used in the curve fits for 
each region (fig. 9). The division lines are located 
at log jo (p/ Po) = —4.5 and log 10 (p/p«) = -0.5. 



Figure 9. Division of curve fit range by density. 


In order to ensure continuity of the dependent vari- 
ables across these two division lines the following 
technique was adopted. If the choice of indepen- 
dent variables yields a point within a specified band 
about either of these division lines, the dependent 
variable is linearly interpolated between the values 
obtained at the endpoints of the band. The coeffi- 
cients for all the curve fits have been tabulated in 
appendix A. In appendix B, a master program which 
handles all the thermodynamic computations is de- 
scribed and a reference is cited for a listing of the 


p = pe{ z ,-\) (‘ 25 ) 

The general form of the equation used for 7 is 

2 

5 = d\ 4* a 2 Y 4 a 3% 4 d 4 Y Z 4 a$Y 
4 a§Z 2 4 ajY 2 Z 4 ag Y Z 1 4 agT 3 
4 aio Z 3 4 (an 4 a\ 2 Y 4 ai3^ 

4 oub Z 4 aisb 2 4 fli6^ 2 + a l7l / Z 
4 a 18 yZ 2 4 a i9 y 3 4 a 20 Z 3 )/[l 4 exp(a 2 i 
4 a22y 4 a 2 3Z 4 Q2 4 Y Z)] (26) 


where Y — log m(p/Po) and Z = log ^{e/RT 0 ). The 
units for p are kg/m 3 and the units for e are m /s~. 
It should be noted that not all the terms appearing in 
the above equation are used over the complete range 
of e and p. 


a = a(e, p) 

An exact expression for the speed of sound a in 
terms of 7 was derived by Barnwell (ref. 13) and may 
be written as 



Since complete bicubic polynomials are used for 
/i(y, Z) and / 2 (F,Z) - f\ [Y,Z) in equation (26) 
for 7, equation (27) is used directly for the corre- 
lation a = a (e,p) without further corrections, un- 
like in references 10 and 15. The expressions for 
ami ( are presented in appendix A. 


T = T(e, p) 

In the calculation of T — T(e,p), the pressure 
is first determined with equation (25), and then the 
temperature is calculated with the equation 

log 10 (T/T 0 ) = hi 4 b 2 Y 4 63Z 4 6 4 yZ 4 b 5 Y 2 4 b 6 Z 2 
4 6 7 y 2 Z 4 6 g yZ 2 4 b 9 Y 3 4 6 10 Z 3 
4 (61 1 46j 2 y + b\sZ + b\ 4 Y Z 
4 6i 5 y 2 4 fcieZ 2 4 b l7 Y 2 Z 4 bi$Y Z 

4 6i 9 y 3 + ^2oZ 3 )/[l4exp(6 21 

4 b 22 Y 4 b 2 sZ 4 b 24 Y Z)] (28) 


where Y = logp) (p/p 0 ), A — log. o(p/Po)i and Z — 
X - y. The units for p are N/m* and the units for 
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7 are K. 1 he coefficients b\ to by \ are determined in 
such a way as to compensate for the errors incurred in 
the initial calculation of pressure with equation (25). 

h = h(p,p) 

For the correlation of b = h(p,p). the ratio 
7 — h/t- U curve fitted as a function of p and p so 
that // can be calculated from 

h = (p/p)[l/(l ~ 0 ] (29) 

l he general form of the equation used for 4 is 

+ c 6 Z 2 + c 7 Y 2 Z + c 8 YZ 2 
4 c 9 y 3 4 ci 0 Z 3 + (cn 4- c\ 2 Y 
+ c i3^ 4 c u Y Z 4 c 1 5 2 -|- C\§Z 2 
+ c 17 Y 2 Z + c 18 YZ 2 +c 19 Y :i 
+ c 20^ 3 )/[1 ± exp(r 2 i 4 C 22 Y 
+ C 23 Z 4- C 24 F Z ) J (30) 

where ) = log 10 (/>//>„), A = log i0 (p/p () ), and Z = 
A - V. For the correlations p = p(e.p) and // = 
//(/j.p). where is the variable curve fitted, an even 
transition function is used to model the transition 
between the perfect-gas equation and the remain- 
der of the curve fit in the lowest density region 
( — 7.0 < ]°&m (p/Po) ^ —4.50). This yields a more 
accurate fit than an ordinary bicubic curve without 
any transitions. 

T = T(pp) 

The general form of the equation used for the 
correlation T — T(p,p) is 

] °S\o( T /To) = d\ 4 - (I2Y 4 d^Z 

4 d 4 YZ 4d! 5 y 2 4 d 6 Z 2 4 d 7 Y 2 Z 
4 ^8^ Z 2 4 dgY ^ 4 d W Z 3 
+ (^11 4d 12 4 + di$Z 4 d\^Y Z 
+ d 15 Y 2 + d l6 Z 2 + d 17 Y 2 Z 

4 d\ 8 Y Z 2 4 digY 3 4 ^20^ 3 )/[l 4 exp(d2l 

4 + ^23-^ + <^24FZ)] (31) 

where V = log i 0 {p/p„), X = log m {p/p„), and 

Z = X-Y. 


s — a(e, p) 

For the correlation of s = .s(e, p), the general form 
of the 1 equat ion used is 

— — i \ 4 <:>) -4 i’:\Z 4- e\Y Z + e^Y 2 4- c^Z 2 
+ ^Y 2 Z -f e s YZ 2 + e 9 V' : * + e w Z :i (32) 

where Y - log w(pjpo) and Z - log j0 (e//?T o ). The 
units for s are nr/s~-K. As is evident from equa- 
tion (32). (irabau transition functions are not neces- 
sary for this curve fit . 

p = pips) 

Unlike the preceding curve fits in which density 
is one of the independent variables, the domain of 
the curve fit p — p(p, .s), as well as the curve fits 
e — e (/4 , s) and a = a(p 4 $), cannot be divided into 
subdomains on the 1 basis of density. For reasons of 
accuracy, it is necessary to subdivide the domain in 
terms of * as shown in figure 10. 



Figure 10. Division of curve fit range by entropy. 


The general form of the equation used for the 
correlat ion of p -- p(p , s) is 

lo Sio {p/po) = f\ 4 f 2 Y 4 / 3 Z 4 f 4 YZ 4 f 5 Y 2 

^ hZ 2 + f 7 Y 2 Z + f 8 YZ 2 

+ /9 V 3 + fioZ 3 + (/1 1 + /12 V 
+- h 3 Z + f 14 YZ + f 15 Y 2 + f 16 Z 2 
5mY 2 z + fisYZ 2 + /i 9 y 3 
/ 20 ^ 3 )/[l + exp(/ 2 i 

+ S22Y + f 23 z + f 24 x + f 25 Y 2 )} (:«) 

where ! = ^ log 10 (.s/ /?), X = log m(p/p„), and 

Z = X — Y . 'Phe units for s are m 2 /s -K. 




e — e(p, s) 

For t hr correlation of c = e(/x .s), the general form 
of the curve fit equation is 

l ° gio ( e / RTo) - g\ 4 - 92 Y + 93 Z + 94 VZ 
+ 9 *>Y 2 + 96% 2 + gjY 2 Z 
+ g%Y Z 2 4 - ggY 3 4 - 9iqZ 3 
+ (ffl 1 + 

+ ffl5V r2 +S 16 ^ 2 + 0l7^ 2 ^ 

+ 018^ Z 2 + 9l9^ 3 ■+• 02O^ 3 )/[1 + exp(ff21 
+ 322^ + 923 Z + 924% + 325 F 2 )] (34) 

where V = \<>%\o(*/ R), X = log l0 (p/p f> ), and 

/ = -V - v . 

a = a(p, «) 

for the correlation of a = n(p, .s), the general form 
of tlu* equation is 

! ° gio ( a / a o ) = + ^2 F + ^ 3 ^ 4 - / i 4 VZ 

4 - fcsK 2 4 - h 6 Z 2 + hyY 2 Z + / i 8 rz 2 
+ / 19F 3 4 - / iio^ 3 4 - (/ii 1 4 -/ 112 ^ 

4 - 4 - /i j 4 K £ 4 " / i j 5 K 2 

+ / i 16 Z 2 4 - / ii 7 V 2 Z 4 - h\ S YZ 2 
+ ^ 1 9 ^ 3 + ^ 20^ 3 )/[l 4 - exp (/ i 2 l 
-4 ^ 22 ^ 4 - h‘ 2 ^Z 4 - / i 24 -^ 4 - ^ 25 ^ 2 )] ( 35 ) 

wIk'it V = l«gio(«//?). A = log 10 (p/p„), and 
^ — A — V. The units of a are m/s. 

Results and Conclusions 

New. simplified curve fits for the thermodynamic 
properties of equilibrium air were constructed with 
the procedures described in the preceding sections. 
Comparisons of t he curve fits p = p(e, />), a = 
a(^p). T = T(e,p), * = .s(e,/i), T - T(p,p), h = 
p = p(p, .s), c = e(p,'.s), and a = a(p, .s) 
with the original NASA RGAS program are shown 
in figures 1 1 to 19. The following procedure was 
employed in making the comparisons for the first four 
curve fits. First, p and p data were supplied as input 
to the NASA RGAS program and e was computed. 
Then, this a and the original p were used to obtain 
P> 7\ and s from the above curve fits. As a result 
of this procedure, log \o{p/p 0 ) is plotted as one of 
the independent variables in figures 11 to 14. The 
same p and p data used above were also employed in 
the comparisons for the curve fits T = T(p, p) and 
h = h(p, p). 

The met hod adopted for t he comparisons of p — 
p(/h' s ) 1 ^ — e(p, .s), and a — a(p, .s) with the NASA 


RGAS program was quite similar to that for the first 
four curve fits. First, /; and p were supplied to the 
NASA RGAS program, which yielded ,s. This $ and 
the original p were used in the above curve fits to 
obtain p . e, and a. 

The above comparisons are presented graphically 
to provide a qualitative overview of the accuracy of 
the curve fits. However, as figures 11 to 19 indi- 
cate, these graphical comparisons are restricted to 
points lying on 1 1 const ant -density lines ranging from 
10 - ' to 10 3 amagats. I 11 order to ensure the validity 
and accuracy of the curve fits across the entire do- 
main, a more comprehensive accuracy test was car- 
ried out. The new curve fits were compared with 
the NASA RGAS program for relative accuracies at 
approximately 22 000 data points. These test points 
were chosen to span the entire density range from 
10 ' to 10 3 amagats and temperatures varying from 
273 K to 25 000 K. The results of these comprehen- 
sive accuracy checks are presented in tables 1 to 9. 
For the curve fits p = p(e,p), a = a(e.p), T = 
T(e, p), T = T(p,p), and h — h(p, p), comparisons 
with the curve fits of reference 15 are also presented 
in the tables. The first column in the tables rep- 
resents the percentage error in the comparison of 
a property generated by the RGAS program and a 
curve fit. The other columns contain the percent- 
age of points in the test data base, generated by a 
curve fit, which are in error by an amount greater 
than that indicated in column I. The accuracies of 
the present curve fits are substantially improved over 
the accuracies of the previous curve fits appearing in 
reference 15. The somewhat higher percentage er- 
rors in the curve fits with p and ,s as independent 
variables can be attributed to the fact that a line of 
constant «s spans the entire density range, sometimes 
necessitating the use of two Grabau-type transition 
functions. Requiring a minimal mismatch across the 
junctions of these transition functions resulted in a 
relative loss of accuracy. However, these latter curve 
fits are well within the accuracy limits required for 
most engineering applications. 

One of the primary objectives of this research 
was to minimize the discontinuities in the depen- 
dent variables across juncture points of the curve fits 
(fig. 20). Comparisons of the dependent variables at 
juncture points of the curve fits for p = p{e,p), a = 
a(e,p), T = T(e, p), T = T(p, /;), and h — h(p, p) 
are presented in tables 10 to 14. These new curve 
fits showed a substantial improvement in continu- 
ity at the juncture points when compared with 
the previous curve fits. For the curve fits where 
and .s were the independent variables, it was 
very difficult to maintain continuity at the junc- 
ture points. This was due to the manner in which 
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Figure 17. Comparison of curve fits for p = p(p, s). 
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Figure 18. Comparison of curve fits for e = e(p.s). 
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very difficult to maintain continuity at the junc- 
ture points. This was due to the manner in which 
the domain was subdivided to obtain the piecewise 
approximating functions. However, discontinuities 
were kept to a minimum, with average mismatches 
of 2.4 percent for p = p(p, .s), 1.2 percent for a — 
a(p, s), and 2.0 percent for e = e(p,s). 



Figure 20. Example curve fit for p = p(e,p). 

A comparison of the relative computer times re- 
quired for the new curve fit subroutines and the 
NASA RGAS program on the National Advanced 
Systems 9160 computer is given in table 15. The 
new subroutine for determining p = p(e,p), a = 
u(e,p), and T — T(e, p) was 2.4 times faster than 
the NASA RGAS subroutine. The previous subrou- 
tine (ref. 15) for the same curve fits was 3.4 times 
faster than the NASA RGAS subroutine. The new 


subroutine for T = T(p, p) was 2.7 times faster 
than the NASA RGAS subroutine, and the previ- 
ous subroutine (ref. 15) was 3.4 times faster. The 
new subroutine for h = h(p,p) was 3.2 times faster 
than the NASA RGAS subroutine, compared with 
the previous subroutine (ref. 15), which was 4.4 times 
faster. The subroutine for s = s (e, /?) was 10.2 times 
faster than the NASA RGAS program. The new 
subroutines for the curve fits p = p(p,s ), e — 
e(p, s), and a = a(p,s) were approximately 10 times 
faster than the NASA RGAS subroutine. It should 
be noted that the NASA RGAS program requires two 
data files for storage of the cubic interpolation coeffi- 
cients. The fact that these data files are now r on disk 
and not tape has significantly speeded up the NASA 
RGAS subroutine. However, the curve fits still pro- 
vide a substantial improvement in computing time, 
being 2.4 to 10.2 times faster than the table-lookup 
technique. 

In conclusion, the new, simplified curve fits for the 
thermodynamic properties of equilibrium air provide 
substantial reductions in computer time and storage 
while maintaining good accuracy. They can be in- 
corporated into computational fluid dynamics com- 
puter codes in a straightforward manner without the 
need for data files. The improved accuracy of the 
new curve fits permits their use in time-dependent 
flow calculations from start-up to the final steady- 
state solution. In addition, the improved continuity 
of these curve fits permits their use in iterative calcu- 
lations. For example, the new curve fit for h = h(p , p) 
can be employed in the iterative procedure required 
to “fit” a bow shock in equilibrium flow. However, 
the discontinuities which still exist in the entropy 
curve fits may cause difficulties when used in an it- 
erative shock calculation. 

NASA Langley Research Center 
Hampton, Virginia 23665-5225 
May 1, 1987 
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Table 1. Accuracy of p — p(c. p) 

Total nuinber of data ])oints = 22 239 

Current, results: Maximum error — 3.93 percent 

log i o(f>/po) = 4.0; log io (e/ RT„) = 0.28 
r = 1.47 x 10 4 K 

Ref. 15 results: Maximum error = 9.00 percent 

•ogio(*>M>) = -4.5; log l0 (e/RT„) = 2.200 
T = 4.50 X 1 0 :i K 



Table 2. Accuracy of a — a(i . p) 

Total number of data points == 22 239 

Current results: Maximum error = 1.18 percent 

log 10 (fi/Po) = -3.0: log 10 (e//?T fJ ) = 0.018 
T = 2.0 X 10 ‘K 

lief. 15 results: Maximum error = 4.94 percent 

log 10 (p/po) = -7.0; logi 0 (e//?7'„) = 0.279 
T = 1.25 x 10 'K 






Table 3. Accuracy of T — T(c. p) 


Total tiurnber of data points = 22 239 
Current results: Maximum error = 4.36 percent 

! <>glO {p/p<>) = -4.0; log l0 (e/RT (> ) = 3.28 
7 = 1.47 x 10 4 K 

H c-f. 15 results: Maxiimitn error = 8.8 percent 

•og|0 {p/po) = -0.625; log io ( e/ RT„) = 3.255 
r = 2.4 x 1 0 1 K 


Error, 

percent 

Current results, 
percent 

Results from 
ref. 15, percent 

0.5 

34.11 

63.82 

1.0 

10.87 

34.74 

2.0 

.58 

9.51 

3.0 

.10 

2.43 

4.0 

.01 

.59 

5.0 

0 

.19 

6.0 

0 

.09 

7.0 

0 

.04 

8.0 

0 

.02 

9.0 

0 

0 

>10.0 

0 

0 


Table 4. Accuracy of ,s = ,s(e, p) 

Tot al number of data points = 21 975 
Current results: Maximum error = 2.51 percent 

login (p/po) = -0.625; log l0 (e/RT„) = 0.657 
7 = 4.89 x 10 2 K 


Error, 

percent 

Current results, 
percent 

0.5 

49.77 ' 

1.0 

15.95 

2.0 

.56 

3.0 

0 

4.0 

0 

5.0 

0 

6.0 

0 

7.0 

0 

8.0 

0 

9.0 

0 

>10.0 

0 




Table 5. Accuracy of T = T(p, /?) 


Total number of data points = 22 239 
Current results: Maxiimim error = 3.9 percent 

*‘>gi() (f>/po) = — 3.25; log 10 (p/p„) - log j = 2.58 
T = 2.4 x 10*K 

Ref. 15 results: Maximum error = 5.71 percent 

l«gio(v>/po) = -0.625; log ,„(/)//>») - log , „(/>//.>„) = 2.44 
r = 2.5 x io'k 


Results from 
ref. 15, percent 
58^82 
28.75 
4.89 
.96 
.16 
.04 

0 

0 

0 

0 

0 


Table 6. Accuracy of h = h(p< p) 

Total number of data ])oints = 22 239 
Current results: Maximum error — 3.44 percent 

login ip/Po) - -7.0; log l0 (p/p») - log w (p/p„) = 2.60 
T = 1.91 x 10 4 K 

Ref. 15 results: Maximum error = 6.56 percent 

log in (/V P<>) = -4.5; log l0 (p/p„) - log, 0 (p/i>„) = 1.01 
T = 2.47 x 10 3 K 

Results from 
ref. 15, percent 
67A5 
40.36 
13.65 
4.78 
1.56 
.46 
.16 

0 
0 
0 
0 




I 




Table 7. Accuracy of p — p(j>. s) 

Total number of data points = 21 030 
Current results: Maximum error = 7.58 percent. 

•oglO (p/po) = -6.625; log 10 (e//?T o ) = 3.30 
T = 1.42 x 10 4 K 



Table 8. Accuracy of e = e(p, s) 

Total number of data points = 21 030 
Current results: Maximum error = 4.5 percent 

•og ID (p/Po) = 2.875; log W (e/RT 0 ) = 2.85 
T = 2.46 x 10 4 K 








Table 10. Comparison of Variables at Juncture Points for p — p(c,p) 


density 

ratio, 

Point A* 

Point B* 

Point C* 

Point D* 

Point E* 

p/po 

Lower 

Upper 

Lower 

Upper 

Lower 

Upper 

Lower 

U pper 

Lower 

Upper 

io- 7 

1.79 x IO- 2 

1.81 x IO -2 

' 7.32 x 10“ 2 

7.38 x 10“ 2 

1.90 x 10“ 1 

1.90 x 10“ 1 

8.72 x 10“ 1 

8.72 x 10“ 1 

2.62 x 10° 

I 2.63 x 10' 

10“ 6 

1.79 x 10“ 1 

1.80 x 10“ 1 

7.81 x 10“ 1 

: 7.78 x 10“ ] 

! 

2.01 x 10° 

2.03 x 10° 

9.53 x 10° 

9.63 x 10° 

2.86 x IO 1 

2.89 x 10 

10“ 5 

1.79 x 10° 

1 .80 x 10° 

1 8.17 x 10° 

8.19 x 10° 

2,16 x 10 1 

2.18 x 10 1 

2 

1.05 x 10 2 

1.06 x 10 2 

3.15 x 10 2 

i 3.16 x 10 

IO" 4 

1.80 x 10 1 

l .81 x 10 1 

8.67 x IO 1 

8.70 x 10 1 

2.40 x 10 2 

2.43 x 10 2 

9.78 x 10 2 

9.79 x IO 2 

1.80 x 10 2 

1.81 x 10 

10“ 3 

1.80 x 10 2 

1.81 x 10 2 

9.12 x 10 2 

9.13 x 10 2 

2.61 x 10 3 

2.63 x IO 3 

1.09 x 10 4 

1.09 x 10 4 



10- 2 

1.80 x IO 3 

1.81 x 10 3 

9.51 x 10 3 

9.51 x 10 3 

2.83 x 10 4 

2.84 x 10 4 

1.23 x 10 5 

1.23 X 10 5 



10 _1 

] .80 x 10 4 

1.81 x 10 4 

9.80 x 10 4 

9.81 x 10 4 

3.08 x 10 5 

3.08 x 10 5 

1.39 x 10 6 

1.39 x 10 6 



10° 

1.80 x 10 s 

1.81 x 10 5 

1.36 x 10 € 

1.36 x 10 6 

4.11 x 10 6 

4.11 x 10 6 





10 1 

1.80 x 10 6 

1.81 x 10 6 

! 

1.41 x 10 7 

1.41 X 10 7 

4.50 x 10 7 

4.54 x 10 7 





10 2 

1 .80 x 1 0 7 

1.81 x 10 7 

1.45 x 10 8 

1.43 x 10 8 

4.91 x 10 8 

5.00 x 10 8 





10 3 

1.80 x 10 s 

1 .83 x 10 8 

1.48 x 10 9 

1.45 x 10 9 

5.42 x 10 9 

5.53 x 10 9 






*See figure 20 for curve breaks. 


Table 11. Comparison of Variables at Juncture Points for a — a(e,p) 


)ensity 

ratio, 

Point A* 

Point B* 

Point C* 

Point D* 

Point K* 

pfpo 

Lower 

Upper 

Lower 

Upper 

Lower 

Upper 

Lower 

U pper 

Lower 

Upper 

10~ 7 

440 

441 

769 

790 

1250 

1260 

2718 

2733 

4731 

4715 

10~ 6 

440 

439 

808 

814 

1291 

1307 

2857 

2871 

4983 

5016 

10“ 5 

440 

438 

831 

841 

1343 

1359 

3021 

3029 

5259 

5287 

10 4 

441 

440 

869 

874 

1429 

1441 

2923 

2925 



10^ 3 

441 

440 

902 

904 

1498 

1506 

3115 

3116 



10“ 2 

441 

440 

932 

932 

1573 

1578 

3337 

3341 



10 1 

441 

441 

957 

957 

1655 

1656 

3596 

3602 



10° 

442 

441 

1120 

1118 

1924 

1924 





10 l 

442 

440 

1149 

1145 

2027 

2039 





10 2 

442 

440 

1171 

1164 

2141 

2166 





10 3 

442 

441 

1188 ! 

1179 

2287 

2312 






See figure 20 for curve breaks. 





Table 12. Comparison of Variables at Juncture Points for T — T(e t p) 


Point B* 


Lower Upper 

486 481 


Lower Upper 

2U2 2091 


Lower Upper 

4033 4034 


Lower Upper 

7 868 7 869 


4283 4284 


2312 2312 


10 364 10 319 

11190 11177 


5307 5326 


2413 2416 


6585 6595 


See figure 20 for curve breaks. 


Table 13. Comparison of Variables at Juncture Points for T 7 T (p,p) 


Density 

ratio, 

Point A* 

Point B* 

Point C* 

Poin 

D* 

p/po 

Lower 

Upper 

Lower 

Upper 

Lower 

Upper 

Lower 

Upper 

Tir 7 

486 

482 

2089 

2089 

4025 

4033 

7 864 

7 838 

10“ 6 

486 

482 

2165 

2165 

4281 

4281 

8 470 

8481 

10“ 5 

486 

484 

2242 

2242 

4549 

4554 

9146 

9146 

10 -4 

486 

482 

2310 

23 L0 

5064 

5042 

10 796 

10 746 

10 -3 

486 

481 

2363 

2363 

5386 

5376 

1 1 793 

11682 

10“ 2 

486 

481 

2404 

2404 

5690 

5701 

12 742 

12 679 

10 _1 

486 

482 

2402 

2402 

5968 

5998 

13671 

13 687 

10° 

486 

482 

2700 

2700 

6248 

6267 



10 1 

486 

482 

2706 

2710 

6585 

6598 



10 2 

486 

483 

2711 

2712 

6950 

6959 



10 3 

486 

483 

2713 

2713 

7309 

7319 




*See figure 20 for curve breaks 



Table 14. Comparison of Variables at Juncture Points for h — h(p,p) 


isity 

tio, 

l^int A* 

Point B* 

Point C* 

Point D* 

Po 

Lower 

Upper 

Lower 

Upper 

Lower 

Upper 

Lower 

Upp 


0.346 x 10* 

0.346 x 10 6 

0.282 x 10 7 

0.285 x 10 7 

0.160 x 10 8 

0.159 x 10 8 

0.997 x 10 8 

0.997 

“ 6 

.346 

.346 

.253 

.254 

,138 

.138 

.890 

.890 

“ 5 

.346 

.346 

.233 

.235 

.120 

.122 

.793 

.792 

-4 

.346 

.346 

345 

.345 

.247 

.247 

.812 

.813 

-3 

.346 

.346 

.314 

.315 

.214 

.214 

.720 

.721 

-2 

.346 

.346 

.296 

.296 

.186 

.186 

.646 

.646 

- 1 

.346 

.346 

.288 

.288 

.164 

.164 

.590 

.591 

0 

.345 

.345 

.386 

.387 

.201 

.202 



1 

.345 

.345 

.377 

.380 

.180 

.181 



2 

.345 

.345 

.374 

.376 

.166 

.166 



3 

.345 

.345 

.374 

.374 

.156 

.156 




See figure 20 for curve breaks. 


Table 15. Comparison of Computer Times 


Computer time, s, for 


Curve fit 

V = P(e, p ) 
a = a(e, p) 

T = T(e,p) 

» = - s (e, p) 
T = T(p, p) 

h = h(p, p) 

P = P{P, ») 
e = e{p, ») 

a = n(r). s) 


Number of 
data points 


10661 

10661 
9921 
9921 
3 038 
3 
3 


Old 

subroutine 
(ref. 15) 


New 

subroutine NASA RGAS 




Appendix A 
Curve Fit Coefficients 

P = P(e.p) 

Tin* coefficients a\, a?- . . u>q and the 1 proper 
sign before the exponential term of the Crabau tran- 
sition in equation (26) are given in tables Al to A3. 
Table Al is for the density range 1 -7.0 < Y < —4.5. 
table A2 is for —4.5 < V < —0.5, and table A3 is for 
—0.5 < V < 3.0. where Y = log'io(p/p»)’ 

The following linear interpolation technique was 
adopted for all the curve fits where density was 
one of the 1 independent variable's. In general, for 
/ = f(YZ), where / is the despondent variable, 
V = log io(p/p<>)' and Z is the second indepen- 
dent variable (either internal energy or pressure), if 
|r - (—4.5)| < 2.5 x KT 2 , then 

f(\\Z) = /( -4.475, Z) + [/(— 4.475, Z) 

~f( -4.525, Z)} 

x (V + 4.525)/0.05 (Al) 

If |V - (-0.5)| < 5.0 x 10' :i . then 

f(Y.Z) = /(— 0.495, Z) + [/(— 0.495, Z) 

- /(- 0.505, Z)\ 

x (V + 0.505 )/().()! (A2) 


a = a(e, p) 

The exact expression for a was given in equa- 
tion (27). The expressions for (757^7) and (777^)^ 
are 1 given below: 

( 


In p / e 


1 4 

In 10 <h 


where 


Oj O n 

— — — ~ a 2 a & ^ z f ‘ 2 or^Y a 2 a r j\ z a ag z 4 3(U) y ~ 

+ (« 1 2 + «i 4 Z 4 2a 15 V -f- 2a 1 j Y Z 4 cl j g Z^ 4- 3ajgF^)/ 
[1 ± exp(a 21 4 a 2 2 F F a 2 3 Z F a 24 V' Z)j 
4 ( a 1 1 F a 1 2 F + aj3Z + ai4yZ + a]5V^4ajf5 Z^ 

+ a 1 7 F 2 Z 4 a^gF Z'* F o 1 9 y ^ 4 a 2 ()Z‘*)(tt 22 + a 24 Z) 
[exp(a 21 4 a 22 V 4 a 23 Z 4 a 24 FZ)/ 

[l±exp(a 21 4 a 22 y 4 a 2 3^ + «24^^)S 2 (A4) 


/ \ 1 Hi 

\() In e ) f) In 10 ()Z 


(A5) 


where 

~ «3 + <mF ) 2«()Z 4 ayT 2 F 2agFZ 4 3aioZ 2 
dZ 

+ (a 13 + a \ A Y 2 a) 6 Z F 2 ai 7 F 2 4 2 ajgFZ 
+ 3 a 20 Z 2 ) / [1 i: exp ( a 2 1 4 a 2 2 F 
+ a 23 Z 4 a Y Z}] (ai 1 + ai 2 F + a^Z 
4 ai 4 Y Z F (Z15V 2 F ajgZ 2 F a]7F 2 Z 
F a 18^ <^ 2 F a 1 9 V 3 F «20^ 3 )( a 23 + a 24T) 

[exp(a 2 1 -*• a. 22 F F «23^ F a 2 4TZ)]/ 

[l±exp(« 21 \ a 22 F F a23^ + q 24TZ ]) 2 (A6) 

The coefficients uj. a 2 , . . . , a 24 are presented in 
tables A 1 to A3. 


r = T(c,p) 

Coefficients /q . /q>, . fr 2 4 are presented in 
tables A l, A5, and A0, where equation (28) gives 
the form of the curve 1 iit. 


h = h{p, p) 

The equation of the curve fit is given by equa- 
tion (30). The coefficients rq, c 2 , . .., c 2 .j and the 
sign before the exponent of the Grabau transition 
function are 1 presented in tables A7, A8, and A9. 

T = T(p, p) 

The coefficients dj, r/ 2 , .... d 2 .j of the curve 
fit, equation (31). an* presented in tables A10. All. 
and A 12. 

s — s(e. p) 

The coefficients rq, e 2 . cqo of the curve fit, 
equation (32). are 1 presented in table A 13. 

P = p{p») 

The general form of the curve fit is given by 
equation (33). The coefficients /j, / 2 . ..., / 2 r, are 
presented in t able All. 

e — e(p. s) 

The coefficients rq, r/ 2 , . . . , g 2 r> of the curve fit, 
equation (34), are presented in table A15. 

a = a(p , 5) 

The curve lit is giv<'n by equation (35). The coef- 
ficients h | , h'>, - - • . /e_>r, are presented in table A 16. 



Table Al. Coefficients for p — p(e , p) and a — a(e,p) for —7.0 < 1 
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02-1 0 0 8.259346/7-01 2.009706E0O -.364 1774 £01 




Tabic A2. Coefficients for p — p(e,p) and a 
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Table A3. Coefficients for p = p(e,p) and a 
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Table A4. Coefficients for T — T(e, p) for —7.0 < Y 
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For Z < 0.25. T = p(e.p)/pR. 




Table A5. Coefficients for T = T(e.p) for -4.5 < Y 
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For Z < 0.25, T — p(e, p)/pR. 




Table A 6. Coefficients for T — T(e,p) for —0.5 < Y < 3.0 
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For Z < 0.25, T — p(e,p)/pR . 




Table A7. Coefficients for h = h(p.p) for 


^ St) 


rr ® 


oooooooooo 

°?o|og^§ 2 


lC ’ — i — < 

o o o 

^ to to 

VI ^ o 
Vl x a 
tsa 10 cc 
N r- O 
\/ -h 


^ 

o o o 


o © o o 

tO tO to to 

co © ic 




SO-^MON^H 

cot^Y ^ 10 — ^ 06 06 

I 1 II i 


t- t- r- o 

CD CO O tF 
GO 00 Tf N 


— < O 

o o o 
tO tq tq 

lO CO O-l 

Tf Ci o 

t - CD <M 
tO N 

o © o 


OCJhc^o^O-MhhM — I <>j ,— i 

ooooooooooooooo 


kO CM CM 

00 o o 
o to Lq 

VI o o 

cO »— i 
iO o 


CM CM o 

o o o 
tq to 143 

© O' CM 
© O CO 
-H ! CO tO 
h» T* 00 
X lO CM 


CO^CMCMCMCMO^O 


ooooooooo 


to tq Gq 

CO 00 CO 
CO N O 
(M C*0 N 


co tp co °c cm 
X O 0 J) - LO 
. OO 

C! X ^ * CM CM 

I I 


tq cq to tq to tj 

^ co lO O X c 

CO CM ^F ^ CO h 

■© cm cm r 


co CO t"- co 

r- x 


o o o o 
o © o © 
to tq tq cq 
o o o o 


N LO M x 
N O Cl N O CO 


o § 

C ^ 
o zo 

V /I X © o o o o o 


cl ci c! ci ci 


33 




Table AS. Coefficients for h = h(p. p) for 
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Table A9. Coefficients for h — h(p, p) for —0.5 < Y 


q 

f 

VI 


N 

V 

o 


H M O ^ 

o o o o 

fcj § b) 

05 N ^ 
C5 rr m C5 

or Str ,-r, 


CM 00 


o ^ 

2 cm 52 

T— ^ . PC' 

CM | i-H 

I 1 


^ 

o o o 

cq cq tq 

CO TP O 
O CM 
O M C5 
O O h- 
^ 05 0© 


O CM 

o o 

cq 


O o 


to S 

S ss 

on 

O H 

05 • 

^ 00 


I I 


Cq c *5 

oo 

O 00 

Tf 

CM 00 
io a 

05 CM 

I 


o 

o 

CM 

co 

co 

r- 

q 

00 

I 


CM CM 

o o 

cq <1 

co 00 

o o 

e- tt 

00 O 'X oc 
CO CM CO 


^ CM 

o o 

&b 9 


o 


CO 


co 

CO 

r- 

I 


N 


s 


VI 

In 

v 

IS 


0 

I 

tl 

CO' 

lC 

r- 

co 

1 cm 

I 


— * CM 

o o 
cq Cq 
co 1 -t 
oo 

CM lC 
CM TJH 

e- q 

cm — 


CM CM 

q o o 

tq tq bq 

i— i lO CM 
^ CO CO 
CM 05 05 
CM -h 00 
q 00 CO 

co 


O O 

o o o 

^ ^ tq 

— o q 

-H CO g 
05 O 
X N g 

lo -oo q 

— rvj CO 


I I 


— CM 

o o 

N eq 

LO CO 
CM O 
*-1 05 

LO T— I 

LO P 

1 7 


cm — tt 

q o o 

iqq c^q [qq 

CO — CM 
CON 
C5 X N 
1C o N 
CO 00 05 


3 

CM 

. — i 

CM 

_ 

o 

o 

o 

o 

o 

fe3 

cq 

cq 

cq 

E*1 

CO- 

ro 

CO 

CO 

CM 

lO 

r- 

05' 

'00 

1—1 

— 

lO 

lC 

e- 

CM 

lO 

oo 

00 

i— t 

CO 

t H 

CM 

CM 

lO 

I - 

CM 

OO 

’O 4 


q 

1 


CO 

1 

1 

CM 


s 


VI 

N 

V 

o 


CM O 

o o 

cq cq 

r- t- 

CM 00 
O O0 
05 CM 

q q 
CM 1-5 

I 


— ! CM 

o o 
cq ^ 

CM Zl 

£ ^ 
-7 


o o 

cq N 

O CM 

■CO o 

N iO 
00 *-• 
tT 00 


— < CM O i— * 

o o o o 

gggg 

SggS 

« 


CM — 

o o 

CjJ r ■>, 

S 30 

q ^ 


o o 

cq cq 

00 LO 
05 CO 
00 00 
05 
. * LO 

I ^ 

1 CM 

I 


+ 


— CM — 


o 

VI 

N 


o 

cq 

o 

tt 


ooooooooooooooooo 


o o o 


cm re — i.o o t- 


o 

o 

o 


oi cc -r lo c i- x 

?-> cj cj, cj G G 


+ 


bp 

tn 


35 




Table A 10. Coefficients for T = T(p.p) for -7.0 < Y 
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I 


For Z < 0.25. T = p/pR. 




Table All. Coefficients for T = T(p.p) for - 4.5 < Y < - 0.5 
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For Z < 0 . 25 . T = p/pR. 



Table A12. Coefficients for T — T(p,p) for -0.5 < 1 
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Table AM. Coefficients for p ~ p(p,s) 
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Table A 15. Coefficients for e 
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Table A 16. Coefficients for 
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Appendix B 
Master Program 

All the curve fits developed in this study have 
been incorporated into a single master program 
called TGAS. This master program permits the user 
to select the desired curve fit from a menu of pos- 
sibilities. The calling statement for this subroutine 
is 

CALL TGAS (P , RHO , E , H , T , A , S , NTGAS) 
where 

P pressure, N/m 2 

RHO density, kg/m'* 

E specific internal energy, m 2 /s 2 

H specific enthalpy, m 2 /s 2 

T temperature, K 

A speed of sound, m/s 

S specific entropy, nr /s-K 

NTGAS integer flag to be set 

by the user for selection 
of the appropriate curve 
fit as fol low's: 


NTGAS = 0: p = p(e, p) 
NTGAS = 1: p = p(e,p), 
a = a(e,p) 
NTGAS = 2: p = p(e, p), 

T = T(e,p) 

NTGAS = 3: p = p(e,p), 

a = fl(e, p), 

T = T(e, p) 

NTGAS = 4: s = ,s(e,p) 
NTGAS = 5: T = T(p , p) 
NTGAS = 6: h . = h(p , p) 
NTGAS = 7: p = p(p, s) 
NTGAS = 8: e = e(p, s) 
NTGAS = 9: a = a{p,s) 


The curve fits for p = p(e, p), a = a(e, p), and T = 
T(e, p) have been placed in the single subroutine 
TGAS1. Subroutine TGAS2 computes s = .s(e,p), sub- 
routine TGAS3 computes T = T(p,p), and subrou- 
tine TGAS4 computes h = h(p,p). The curve fits for 
p = p(p, ,s), e = e(p, s , ), and a -= a{p,s) have been 
placed in subroutines TGAS5, TGAS6, and TGAS7, re- 
spectively. The subroutines TGAS1 to TGAS7 can be 
used in a “stand-alone” manner if so desired, inde- 
pendent of the master program. A FORTRAN listing 
of each subroutine is given in reference 18. 
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